﻿#include <xy/alg/sampler.h>

namespace xy
{

namespace alg
{

using namespace xy::geom;
using namespace xy::math;

std::vector<vec2> poisson_disk_sampler::perform(double uFirst, double uLast, double vFirst, double vLast) const
{
    static std::mt19937 generator(std::random_device{}());
    double U1 = uFirst, U2 = uLast, V1 = vFirst, V2 = vLast;
    double dist = std::max(U2 - U1, V2 - V1) / 100;

    std::uniform_real_distribution<double> uDistribution(U1, U2);
    std::uniform_real_distribution<double> vDistribution(V1, V2);

    double width = std::abs(U2 - U1);
    double height = std::abs(V2 - V1);
    double tol = std::min(width, height) / 100;
    double cellSize = dist / sqrt(2);

    // 网格中的元素可以为空
    std::vector<std::vector<vec2 *>> grid(height / cellSize, std::vector<vec2 *>(width / cellSize, nullptr));

    // 网格索引映射
    auto map2grid = [&](vec2 P, int &i, int &j) {
        std::tie(i, j) = std::make_tuple((P.y() - V1) / cellSize, (P.x() - U1) / cellSize);
        i = std::min(std::max(i, 0), (int)grid.size() - 1);
        j = std::min(std::max(j, 0), (int)grid.front().size() - 1);
    };

    // 第一个随机点
    vec2 firstPoint{uDistribution(generator), vDistribution(generator)};

    std::queue<vec2> processQueue;
    processQueue.push(firstPoint);

    std::vector<vec2> samplePoints;
    samplePoints.push_back(firstPoint);

    int i, j;
    map2grid(firstPoint, i, j);
    grid[i][j] = &firstPoint;

    // 生成到 P 距离在 1-2 倍 dist 之间的随机点
    auto genPointAround = [&](vec2 P, double dist) -> vec2 {
        static std::uniform_real_distribution<double> r(0, 1);
        double r1 = r(generator), r2 = r(generator);
        double radius = dist * (r1 + 1);
        double angle = 2 * numeric::pi * r2;
        return vec2{P.x() + radius * std::cos(angle), P.y() + radius * std::sin(angle)};
    };

    // 判断点是否在参数域中
    auto inRect = [&](vec2 P) -> bool { return !(P.x() < U1 || P.x() > U2 || P.y() < V1 || P.y() > V2); };

    // 判断采样点所在单元格的 4 个角点为中心，半径为 dist 的圆内是否有其它点
    // 简单起见，判断周围 25 个格点
    auto inNeighbourhood = [&](vec2 P, double dist) -> bool {
        int p, q;
        map2grid(P, p, q);

        for (int i = -2; i < 3; i++)
        {
            for (int j = -2; j < 3; j++)
            {
                int ii = std::min(std::max(p + i, 0), (int)grid.size() - 1);
                int jj = std::min(std::max(q + j, 0), (int)grid.front().size() - 1);
                if (grid[ii][jj] != nullptr && grid[ii][jj]->distance(P) < dist)
                    return true;
            }
        }
        return false;
    };

    while (!processQueue.empty())
    {
        auto point = processQueue.front();
        processQueue.pop();

        for (int i = 0; i < newPointCount; i++)
        {
            // 第一个点随机生成，增加随机核心数
            vec2 newPoint;
            if (i == 0)
                newPoint = vec2{uDistribution(generator), vDistribution(generator)};
            else
                newPoint = genPointAround(point, dist * std::max(dense(point), tol));

            if (inRect(newPoint) && !inNeighbourhood(newPoint, dist * std::max(dense(newPoint), tol)))
            {
                int p, q;
                map2grid(newPoint, p, q);
                if (grid[p][q] == nullptr)
                {
                    grid[p][q] = &newPoint;
                    processQueue.push(newPoint);
                    samplePoints.push_back(newPoint);
                }
            }
        }
    }
    return samplePoints;
}

} // namespace alg

} // namespace xy
